#simulation function
function Simulate(prim::Primitives, est::Estimands, res::Results, utilities::Vector{Matrix{Float64}}, nsim::Int64; 
    rec::Int64 = 0, cfact::Int64 = 0, dlim::Int64 = 0, decomp::Int64 = 0, rec_pol::Int64 = 0)
    @unpack r, η, β, dbar, cfloor, θ_grid, nθ, a_grid, na, t_grid, nt, e_grid, ne, x_grid, nx, I_grid, W_grid, r_d, d_grid, nd, prob_forb, p_min, p_max, thresh_sub, dbar_baseline, dbar_nolim, dbar_rec = prim
    @unpack γ_grid, λ_0, λ_1, λ_2, λ_3, λ_4, ξ, σ_l, π_0, π_1, π_2, π_3, π_4, π_5, π_6, ω, α, ψ_0, ψ_1, ψ_2, σ_ζ, ψ_3, ψ_4 = est = est

    ###setup
    nrow = nsim * 20 * nt #number of simulations per ventile, one line of data for each age
    ncol = 25 #or whatever
    data_simul = SharedArray{Float64}(nrow, ncol)

    ##reset simulated data, just in case!
    for r = 1:nrow, c = 1:ncol
        data_simul[r, c] = 0.0
    end

    educ_costs, grad_probs, grants_func, grants_merit_func = utilities[1], utilities[2], utilities[3], utilities[4]
    parent_dist_inc_abil, parent_house_dist = utilities[5], utilities[6]
    d_opts, rec_h_shocks = utilities[9], utilities[10]

    #debt forgiveness counterfactual
    d_reduce = 0.0


    #adjust borrowing limit if in recession: 31,000 in 2008 dollars
    if rec_pol == 1 && dlim == 0 && (decomp == 0 || decomp == 3)
        dbar = dbar_rec    
    end

    #removal of borrowing constraints
    if dlim !=0
        dbar = d_opts[dlim, 1]     
    end

    #set up some useful interpolations  
    interp_a = interpolate(a_grid, BSpline(Linear()))  
    interp_d = interpolate(d_grid, BSpline(Linear()))
    interp_p_T = interpolate(res.pol_func_p_τ, (BSpline(Linear()), BSpline(Linear()), BSpline(Linear()))) #total parent transfer function
    
    interp_val_emp = interpolate(res.val_func_emp, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear())))
    interp_val_unemp = interpolate(res.val_func_unemp, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear()), NoInterp()))    
    interp_val_nilf = interpolate(res.val_func_nilf, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear()), NoInterp()))    
    
    interp_pol_emp = interpolate(res.pol_func_emp, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear())))
    interp_pol_unemp = interpolate(res.pol_func_unemp, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear()), NoInterp()))    
    interp_pol_nilf = interpolate(res.pol_func_nilf, (BSpline(Linear()), NoInterp(), NoInterp(), BSpline(Linear()), NoInterp(), BSpline(Linear()), NoInterp()))    

    ####main loop over parent ventiles
    @sync @distributed for i_v = 1:20
        Random.seed!(i_v*10)
        I = parent_dist_inc_abil[i_v, 1]
        I_cfact = I * (1 - .1374) #counterfactual parent income for Hilger elasticity
        
        #distributions of child ability and wealth conditional on ventile
        θ_dist = Normal(parent_dist_inc_abil[i_v, 2], parent_dist_inc_abil[i_v, 3]) #child ability distribution
        dist_H = LogNormal(parent_house_dist[i_v, 2], parent_house_dist[i_v, 3]) #conditional distribution
        dist_Wnh = LogNormal(parent_house_dist[i_v, 4], parent_house_dist[i_v, 5])
        prob_house = parent_house_dist[i_v, 1]

        ###loop over simulations
        for i_sim = 1:nsim       
            #draw child ability and truncate
            θ = rand(θ_dist)
            θ = max(θ, 0.0)
            θ = min(θ, 1.0)
            θ_ind = get_index(θ, θ_grid)

            #re-initialize parent income and introduce layoffs
            I = parent_dist_inc_abil[i_v, 1]
            layoff_draw = rand()
            layoff_probs = utilities[8]  
            if rec == 1 && layoff_draw<layoff_probs[i_v, 4] && (decomp == 0 || decomp == 2)  #if in recession and get a bad layoff draw
                I *= layoff_probs[i_v, 2]
            end
        
            educ_costs_temp = copy(educ_costs)
            gov_exp = 0 #initialize government expenditure on individual
   
            #draw house equity and non-house wealth
            H, H_cfact = 0.0, 0.0 #preallocate housing wealth to zero
            Wnh = rand(dist_Wnh) #draw non-housing wealth

            #draw
            house_draw = rand() #governs whether we have a house         
            if house_draw<prob_house #kid's parents have a house!
                H = rand(dist_H)
                H_cfact = H + 0.25 #normalization by 40k --> 0.25 = 10k
            end

            if rec>0 && (decomp == 0 || decomp == 2) #counterfactual recession
                H *= rec_h_shocks[i_v] #slash house values
            end
            W = Wnh + H #total wealth

            #modify tuition prices in recessions
            if rec == 1 && (decomp == 0)
                educ_costs_temp[2] *= 1.00
                educ_costs_temp[3] *= 1.08
                educ_costs_temp[4] *= 1.03
                educ_costs_temp[5] *= 1.03
            end

            #finanical aid
            grants_aid = zeros(5)
            grants_aid_cfact = zeros(5)
            grants_fed = zeros(5)
            
            for s = 1:4
                grants_aid[s+1] = Grants(I, s, grants_func, W)[1]
                grants_aid_cfact[s+1] = Grants(I_cfact, s, grants_func, W)[1]
                grants_fed[s+1] = Grants(I, s, grants_func, W)[2]
                
                #update grant aid if we're running a counterfactual
                if cfact>0
                    grants_aid[s+1] = Grants(I, s, grants_func, W; cfact = cfact)[1]
                    grants_aid_cfact[s+1] = Grants(I_cfact, s, grants_func, W; cfact = cfact)[1]
                    grants_fed[s+1] = Grants(I, s, grants_func, W; cfact = cfact)[2]
                end
            end
 
            #merit grants
            grants_merit = zeros(5)
            for s = 1:4
                grants_merit[s+1] = Grants_merit(θ, s, grants_merit_func)
            end

            #update pell receipts in recession
            if I<0.820 && rec_pol == 1 && (decomp == 0 || decomp == 3)
                #updating policy
                for s = 1:4
                    temp = grants_aid[s+1] - grants_fed[s+1] + (grants_fed[s+1]*1.25)
                    temp2 = grants_fed[s+1] * 1.25
                    grants_aid[s+1] = temp
                    grants_fed[s+1] = temp2
                end
            end

            #total grants
            grants_total = grants_aid .+ grants_merit
            grants_total_cfact = grants_aid_cfact .+ grants_merit

            #indices
            I_index = get_index(I, I_grid)
            I_cfact_index = get_index(I_cfact, I_grid)
            W_index = get_index(Wnh + H, W_grid)
            W_cfact_index = get_index(Wnh + H_cfact, W_grid)

            #parental transfers
            τ = interp_p_T(I_index, W_index, θ_ind)
            τ_cfact_I = interp_p_T(I_cfact_index, W_index, θ_ind)
            τ_cfact_W = interp_p_T(I_index, W_cfact_index, θ_ind)

            ###compute value of kid schooling options
            ns = length(educ_costs_temp) #number of schooling options
            s_vals, s_vals_cfact_I, s_vals_cfact_W = zeros(ns), zeros(ns), zeros(ns)
             
            prob_emp_hs = 𝚽(π_0)
            prob_emp_drop = 𝚽(π_0)
            prob_emp_sc = 𝚽(π_0 + π_2)
            prob_emp_coll = 𝚽(π_0 + π_3)
            rec_emp_shocks = utilities[7]            

            #adjust initial employment probabilities if in recession and parents paying attention
            if rec == 1 && (decomp == 0 || decomp == 1)
                #prob_emp_hs = 𝚽(π_0 * rec_emp_shocks[1, 1])
                #prob_emp_drop = 𝚽(π_0 * rec_emp_shocks[2, 1])
                prob_emp_hs = 0.7367
                prob_emp_drop = 0.7491
                prob_emp_sc = 𝚽(π_0 * rec_emp_shocks[3, 1] + π_2 * rec_emp_shocks[3, 2])
                prob_emp_coll = 𝚽(π_0 * rec_emp_shocks[5, 1] + π_3 * rec_emp_shocks[5, 3])
            end

            #value of high school option
            a_start_hs = 0.0 #you're on your own!
            d_start_hs = 0.0 #no debt, though!
            a_hs_ind = get_index(a_start_hs, a_grid)
            d_hs_ind = get_index(d_start_hs, d_grid)

            val_work = prob_emp_hs * interp_val_emp(θ_ind, 1, 1, a_hs_ind, 1, d_hs_ind) + (1-prob_emp_hs) * interp_val_unemp(θ_ind, 1, 1, a_hs_ind, 1, d_hs_ind, 1)
            val_nwork = interp_val_nilf(θ_ind, 1, 1, a_hs_ind, 1, d_hs_ind, 1) 

            #working threshold rule (no switching costs!) and expected max utility
            ε_l_thresh = val_nwork - val_work
            val_hs = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
            s_vals[1] = val_hs
            s_vals_cfact_I[1] = val_hs
            s_vals_cfact_W[1] = val_hs

            #loop over other schooling options
            for i_s = 2:ns
                prob_drop, prob_sc, prob_coll = Grad_prob(θ, i_s-1, grad_probs) #probability of graduation outcomes

                ###starting assets/debt for each possible educational outcome. Bounded by zero in either case
                #dropout
                a_start_drop = max((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0)
                d_start_drop = min((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0) * -1 

                #some college
                a_start_sc = max((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0)
                d_start_sc = min((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0) * -1 

                #college
                a_start_coll = max((τ - (educ_costs_temp[i_s] - grants_total[i_s])), 0)
                d_start_coll = min((τ - (educ_costs_temp[i_s] - grants_total[i_s])), 0) * -1 

                #sub/unsubsidized loan check. If parent income is too high, interest accumulates while in college.   
                if I>thresh_sub
                    d_start_drop *=(1+r_d)
                    d_start_sc *= (1+r_d)^2
                    d_start_coll *= (1+r_d)^4
                end

                #debt forgiveness counterfactual
                if cfact == 1
                    d_start_drop = max(d_start_drop - d_reduce, 0)
                    d_start_sc = max(d_start_sc - d_reduce, 0)
                    d_start_coll = max(d_start_coll - d_reduce, 0)
                end

                #indices
                a_d_i = get_index(a_start_drop, a_grid)
                a_s_i = get_index(a_start_sc, a_grid)
                a_c_i = get_index(a_start_coll, a_grid)
                d_d_i = get_index(d_start_drop, d_grid)
                d_s_i = get_index(d_start_sc, d_grid)
                d_c_i = get_index(d_start_coll, d_grid)

                ###starting assets for each possible educational outcome, I cfact###
                #assets
                a_start_drop_cfact_I = max((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s]))/4, 0)
                a_start_sc_cfact_I = max((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s]))/2, 0)
                a_start_coll_cfact_I = max((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s])), 0)
                
                #debt
                d_start_drop_cfact_I = min((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s]))/4, 0) * -1
                d_start_sc_cfact_I = min((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s]))/2, 0) * -1
                d_start_coll_cfact_I = min((τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s])), 0) * -1

                if I_cfact>thresh_sub
                    d_start_drop_cfact_I *=(1+r_d)
                    d_start_sc_cfact_I *= (1+r_d)^2
                    d_start_coll_cfact_I *= (1+r_d)^4
                end

                #indices
                a_d_i_cf_I = get_index(a_start_drop_cfact_I, a_grid)
                a_s_i_cf_I = get_index(a_start_sc_cfact_I, a_grid)
                a_c_i_cf_I = get_index(a_start_coll_cfact_I, a_grid)
                d_d_i_cf_I = get_index(d_start_drop_cfact_I, d_grid)
                d_s_i_cf_I = get_index(d_start_sc_cfact_I, d_grid)
                d_c_i_cf_I = get_index(d_start_coll_cfact_I, d_grid)

                ###starting assets for each possible educational outcome, W cfact###
                #assets
                a_start_drop_cfact_W = max((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0)
                a_start_sc_cfact_W = max((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0)
                a_start_coll_cfact_W = max((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s])), 0)
                
                #debt
                d_start_drop_cfact_W = min((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0) * -1
                d_start_sc_cfact_W = min((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0) * -1
                d_start_coll_cfact_W = min((τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s])), 0) * -1
                
                if I>thresh_sub
                    d_start_drop_cfact_W *=(1+r_d)
                    d_start_sc_cfact_W *= (1+r_d)^2
                    d_start_coll_cfact_W *= (1+r_d)^4
                end

                #indices
                a_d_i_cf_W = get_index(a_start_drop_cfact_W, a_grid)
                a_s_i_cf_W = get_index(a_start_sc_cfact_W, a_grid)
                a_c_i_cf_W = get_index(a_start_coll_cfact_W, a_grid)
                d_d_i_cf_W = get_index(d_start_drop_cfact_W, d_grid)
                d_s_i_cf_W = get_index(d_start_sc_cfact_W, d_grid)
                d_c_i_cf_W = get_index(d_start_coll_cfact_W, d_grid)

                #######Baseline#######
                ###utility from dropout: threshold, then compute
                val_work = prob_emp_drop * interp_val_emp(θ_ind, 1, 1, a_d_i, 2, d_d_i) + (1-prob_emp_drop) * 
                (interp_val_unemp(θ_ind, 1, 1, a_d_i, 2, d_d_i, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 1, 1, a_d_i, 2, d_d_i, 1) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 1, 1, a_d_i, 2, d_d_i, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 1, 1, a_d_i, 2, d_d_i, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_drop = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals[i_s] += β * val_drop * prob_drop

                #utility from some college
                val_work = prob_emp_sc * interp_val_emp(θ_ind, 2, 1, a_s_i, 3, d_s_i) + (1-prob_emp_sc) * 
                (interp_val_unemp(θ_ind, 2, 1, a_s_i, 3, d_s_i, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 2, 1, a_s_i, 3, d_s_i, 2) * prob_forb) 
                val_nwork = interp_val_nilf(θ_ind, 2, 1, a_s_i, 3, d_s_i, 1) *(1-prob_forb) + interp_val_nilf(θ_ind, 2, 1, a_s_i, 3, d_s_i, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_sc = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals[i_s] += β^2 * val_sc * prob_sc

                #utilty from college
                val_work = prob_emp_coll * interp_val_emp(θ_ind, 3, 1, a_c_i, 5, d_c_i) + (1-prob_emp_coll) * 
                (interp_val_unemp(θ_ind, 3, 1, a_c_i, 5, d_c_i, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 3, 1, a_c_i, 5, d_c_i, 2) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 3, 1, a_c_i, 5, d_c_i, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 3, 1, a_c_i, 5, d_c_i, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_coll = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals[i_s] += β^4 * val_coll * prob_coll

                #######Counterfactual Income#######
                #dropout
                val_work = prob_emp_drop * interp_val_emp(θ_ind, 1, 1, a_d_i_cf_I, 2, d_d_i_cf_I) + (1-prob_emp_drop) * 
                (interp_val_unemp(θ_ind, 1, 1, a_d_i_cf_I, 2, d_d_i_cf_I, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 1, 1, a_d_i_cf_I, 2, d_d_i_cf_I, 1) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 1, 1, a_d_i_cf_I, 2, d_d_i_cf_I, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 1, 1, a_d_i_cf_I, 2, d_d_i_cf_I, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_drop = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_I[i_s] += β * val_drop * prob_drop

                #utility from some college
                val_work = prob_emp_sc * interp_val_emp(θ_ind, 2, 1, a_s_i_cf_I, 3, d_s_i_cf_I) + (1-prob_emp_sc) * 
                (interp_val_unemp(θ_ind, 2, 1, a_s_i_cf_I, 3, d_s_i_cf_I, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 2, 1, a_s_i_cf_I, 3, d_s_i_cf_I, 2) * prob_forb) 
                val_nwork = interp_val_nilf(θ_ind, 2, 1, a_s_i_cf_I, 3, d_s_i_cf_I, 1) *(1-prob_forb) + interp_val_nilf(θ_ind, 2, 1, a_s_i_cf_I, 3, d_s_i_cf_I, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_sc = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_I[i_s] += β^2 * val_sc * prob_sc

                #utilty from college
                val_work = prob_emp_coll * interp_val_emp(θ_ind, 3, 1, a_c_i_cf_I, 5, d_c_i_cf_I) + (1-prob_emp_coll) * 
                (interp_val_unemp(θ_ind, 3, 1, a_c_i_cf_I, 5, d_c_i_cf_I, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 3, 1, a_c_i_cf_I, 5, d_c_i_cf_I, 2) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 3, 1, a_c_i_cf_I, 5, d_c_i_cf_I, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 3, 1, a_c_i_cf_I, 5, d_c_i_cf_I, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_coll = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_I[i_s] += β^4 * val_coll * prob_coll

                #######Counterfactual Wealth#######
                #dropout
                val_work = prob_emp_drop * interp_val_emp(θ_ind, 1, 1, a_d_i_cf_W, 2, d_d_i_cf_W) + (1-prob_emp_drop) * 
                (interp_val_unemp(θ_ind, 1, 1, a_d_i_cf_W, 2, d_d_i_cf_W, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 1, 1, a_d_i_cf_W, 2, d_d_i_cf_W, 1) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 1, 1, a_d_i_cf_W, 2, d_d_i_cf_W, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 1, 1, a_d_i_cf_W, 2, d_d_i_cf_W, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_drop = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_W[i_s] += β * val_drop * prob_drop

                #utility from some college
                val_work = prob_emp_sc * interp_val_emp(θ_ind, 2, 1, a_s_i_cf_W, 3, d_s_i_cf_W) + (1-prob_emp_sc) * 
                (interp_val_unemp(θ_ind, 2, 1, a_s_i_cf_W, 3, d_s_i_cf_W, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 2, 1, a_s_i_cf_W, 3, d_s_i_cf_W, 2) * prob_forb) 
                val_nwork = interp_val_nilf(θ_ind, 2, 1, a_s_i_cf_W, 3, d_s_i_cf_W, 1) *(1-prob_forb) + interp_val_nilf(θ_ind, 2, 1, a_s_i_cf_W, 3, d_s_i_cf_W, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_sc = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_W[i_s] += β^2 * val_sc * prob_sc

                #utilty from college
                val_work = prob_emp_coll * interp_val_emp(θ_ind, 3, 1, a_c_i_cf_W, 5, d_c_i_cf_W) + (1-prob_emp_coll) * 
                (interp_val_unemp(θ_ind, 3, 1, a_c_i_cf_W, 5, d_c_i_cf_W, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, 3, 1, a_c_i_cf_W, 5, d_c_i_cf_W, 2) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, 3, 1, a_c_i_cf_W, 5, d_c_i_cf_W, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, 3, 1, a_c_i_cf_W, 5, d_c_i_cf_W, 2) * prob_forb
                ε_l_thresh = val_nwork - val_work
                val_coll = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                s_vals_cfact_W[i_s] += β^4 * val_coll * prob_coll

                #psychic cost of enrollment
                d_exp = prob_drop * d_start_drop + prob_sc * d_start_sc + prob_coll * d_start_coll    
                ψ_total = ψ_0 + ψ_1 * θ + ψ_2 * θ * (i_s==3 || i_s==5) + ψ_3 * (d_exp>0) + ψ_4 * d_exp

                s_vals[i_s] += ψ_total
                s_vals_cfact_I[i_s] += ψ_total
                s_vals_cfact_W[i_s] += ψ_total

                #check that we don't go over the borrowing constraint
                debt_raw_max = τ - (educ_costs_temp[i_s] - grants_total[i_s])
                debt_raw_max_cfact_I = τ_cfact_I - (educ_costs_temp[i_s] - grants_total_cfact[i_s])
                debt_raw_max_cfact_W = τ_cfact_W - (educ_costs_temp[i_s] - grants_total[i_s])
                
                if (debt_raw_max + d_reduce)<(-dbar) #over borrowing limit. Can't go over limit with forgiven debt
                    s_vals[i_s] = -Inf #not allowed to exceed borrowing constraint!
                end    

                #check that we don't go over the borrowing constraint
                if (debt_raw_max_cfact_I + d_reduce)<(-dbar) #over borrowing limit. Can't go over limit with forgiven debt
                    s_vals_cfact_I[i_s] = -Inf #not allowed to exceed borrowing constraint!
                end    

                #check that we don't go over the borrowing constraint
                if (debt_raw_max_cfact_W + d_reduce)<(-dbar) #over borrowing limit. Can't go over limit with forgiven debt
                    s_vals_cfact_W[i_s] = -Inf #not allowed to exceed borrowing constraint!
                end    
            end

            #enrollment probabiltiies
            s_vals_exp = exp.((1/σ_ζ).*s_vals)
            s_vals_exp_cfact_I = exp.((1/σ_ζ).*s_vals_cfact_I)
            s_vals_exp_cfact_W = exp.((1/σ_ζ).*s_vals_cfact_W)

            #cfact probability of enrolling
            prob_coll_baseline, prob_coll_cfact_I, prob_coll_cfact_W = 0, 0, 0     
            prob_coll_baseline = 1 - s_vals_exp[1] / sum(s_vals_exp)
            prob_coll_cfact_W = 1 - s_vals_exp_cfact_W[1] / sum(s_vals_exp_cfact_W)
            prob_coll_cfact_I = 1 - s_vals_exp_cfact_I[1] / sum(s_vals_exp_cfact_I)

            #Hilger moment is tricky -- it's probably of enrollment in age 18-22 span, pooled. So need to do some summing.
            prob_enroll_age = zeros(5)
            prob_enroll_age_cfact = zeros(5)
            for i_s = 2:5

                prob_drop, prob_sc, prob_coll = Grad_prob(θ, i_s-1, grad_probs) #probability of graduation outcomes

                #probability of schooling option, and probabiliyt of achievemnt outcomes if selected
                prob_s = s_vals_exp[i_s] / sum(s_vals_exp)

                #add to enrollment age likelihoods
                prob_enroll_age[1] += prob_drop * prob_s
                prob_enroll_age[2] += prob_drop * prob_s
                prob_enroll_age[3] += prob_sc * prob_s 
                prob_enroll_age[4] += prob_coll * prob_s
                prob_enroll_age[5] += prob_coll * prob_s

                #now in counterfactual scenario
                prob_s_cfact = s_vals_exp_cfact_I[i_s] / sum(s_vals_exp_cfact_I)
                prob_enroll_age_cfact[1] += prob_drop * prob_s_cfact
                prob_enroll_age_cfact[2] += prob_drop * prob_s_cfact
                prob_enroll_age_cfact[3] += prob_sc * prob_s_cfact 
                prob_enroll_age_cfact[4] += prob_coll * prob_s_cfact
                prob_enroll_age_cfact[5] += prob_coll * prob_s_cfact
            end

            #treatment effects
            prob_diff_cfact_I = (mean(prob_enroll_age) - mean(prob_enroll_age_cfact)) * 100 #with the pooling. I think this will work . . .
            prob_diff_cfact_W = (prob_coll_cfact_W - prob_coll_baseline) * 100

            #draw schooling decisions, outcomes, and shocks
            draws_gumbel = rand(Gumbel(0, σ_ζ), ns)
            s_choice_val = maximum(s_vals.+draws_gumbel) #maximizing value
            s_choice = findmin(abs.(s_vals.+draws_gumbel.-s_choice_val))[2] #index of choice [1-5]
            
            ###now we have our educ choice! Figure out starting indices
            i_e, i_x, i_a, t_start, i_d = 1, 1, a_hs_ind, 1, 1  #this corresponds to the high school case
            
            #now update if the agent went to some sort of college!
            grad_draw = rand() #
            if s_choice!=1
                prob_drop, prob_sc, prob_coll = Grad_prob(θ, s_choice-1, grad_probs) #probability of graduation outcomes

                if grad_draw<prob_drop #dropout
                    t_start = 2
                    a_start_drop = max((τ - (educ_costs_temp[s_choice] - grants_total[s_choice]))/4, 0)
                    d_start_drop = min((τ - (educ_costs_temp[s_choice] - grants_total[s_choice]))/4, 0) * -1 

                    #subsidized loan eligibility check
                    if I>thresh_sub
                        d_start_drop *=(1+r_d)
                    end
    
                    #indices
                    i_a = get_index(a_start_drop, a_grid)
                    i_d = get_index(d_start_drop, d_grid)

                    #govt expenditure
                    gov_exp += grants_fed[s_choice]/4
                    gov_exp += min(d_reduce, max((educ_costs_temp[s_choice] - grants_fed[s_choice] - τ)/4, 0.0))

                elseif grad_draw<prob_drop + prob_sc && grad_draw>prob_drop #some college
                    i_e, t_start = 2, 3
                    a_start_sc = max((τ - (educ_costs_temp[s_choice] - grants_total[s_choice]))/2, 0)
                    d_start_sc = min((τ - (educ_costs_temp[s_choice] - grants_total[s_choice]))/2, 0) * -1 

                    #subsidized loan eligibility check
                    if I>thresh_sub
                        d_start_sc *= (1+r_d)^2
                    end
                
                    i_a = get_index(a_start_sc, a_grid)
                    i_d = get_index(d_start_sc, d_grid)
                    
                    #govt expenditure
                    gov_exp += grants_fed[s_choice]/2
                    gov_exp += min(d_reduce, max((educ_costs_temp[s_choice] - grants_fed[s_choice] - τ)/2, 0.0))

                elseif grad_draw>prob_drop + prob_sc #bachelor's
                    i_e, t_start = 3, 5
                    a_start_coll = max((τ - (educ_costs_temp[s_choice] - grants_total[s_choice])), 0)
                    d_start_coll = min((τ - (educ_costs_temp[s_choice] - grants_total[s_choice])), 0) * -1

                    if I>thresh_sub
                        d_start_coll *= (1+r_d)^4
                    end

                    i_a = get_index(a_start_coll, a_grid)
                    i_d = get_index(d_start_coll, d_grid)

                    #govt expenditure
                    gov_exp += grants_fed[s_choice]
                    gov_exp += min(d_reduce, max((educ_costs_temp[s_choice] - grants_fed[s_choice] - τ), 0.0))
                end
            end
            
            emp_shocks = rand(nt) #preallocate vector of employment draw shcoks
            forb_shocks = rand(nt) #preallocate vector of forebearance draw shocks
            l_shocks = rand(Normal(0, σ_l), nt) #labor supply preference shocks            
            p = 0 #previous participation state, for keeping track of switching costs

            #####Simulate lifetime earnings path
            for i_t = t_start:nt #from 18/19/20/22 to 65
                
                #unpack state space
                e, x, a, t, d = e_grid[i_e], x_grid[i_x], interp_a(i_a), t_grid[i_t], interp_d(i_d)
                
                rec_emp_shocks = utilities[7]
                π_0_r = π_0
                π_2_r = π_2
                π_3_r = π_3
                if rec==1 && i_t<6 && (decomp == 0 || decomp == 1)
                    π_0_r = π_0 * rec_emp_shocks[i_t, 1]
                    π_2_r = π_2 * rec_emp_shocks[i_t, 2]
                    π_3_r = π_3 * rec_emp_shocks[i_t, 3]
                end

                ###start of period: figure out labor supply decision
                prob_emp = 𝚽(π_0_r + π_1 * (p==2) + π_2_r * (e==2) + π_3_r * (e==3) + π_4 * x + π_5 * x * (e==2) + π_6 * x * (e==3))

                #hard-coding employment probabilities for high schoolers at start of recession
                if rec == 1 && i_t == 1 && s_choice == 1 && (decomp == 0 || decomp == 1)
                    prob_emp = 0.7367

                end

                #same, but for dropouts
                if rec == 1 && i_t == 2 && s_choice!=1 && i_e == 1 && (decomp == 0 || decomp == 1)
                    prob_emp = 0.7491
                end

                val_work = prob_emp * interp_val_emp(θ_ind, i_e, i_x, i_a, i_t, i_d) + (1-prob_emp) * 
                (interp_val_unemp(θ_ind, i_e, i_x, i_a, i_t, i_d, 1) * (1-prob_forb) + interp_val_unemp(θ_ind, i_e, i_x, i_a, i_t, i_d, 2) * prob_forb)
                val_nwork = interp_val_nilf(θ_ind, i_e, i_x, i_a, i_t, i_d, 1) * (1-prob_forb) + interp_val_nilf(θ_ind, i_e, i_x, i_a, i_t, i_d, 1) * prob_forb

                #baseline case
                ε_l_thresh = val_nwork - val_work

                #update if not in intiial period
                if i_t>t_start && p == 0
                    ε_l_thresh += ξ
                elseif i_t>t_start && p!=0
                    ε_l_thresh -= ξ
                end

                ###determine LFP and employment
                l, m, forb = 0, 0, 0
                if l_shocks[i_t]>ε_l_thresh #chooses to work!
                    l = 1
                    if emp_shocks[i_t]<prob_emp #employed!
                        m = 1
                    end
                end

                #determine forebearance: need to be not employed, have positive debt, and correct shock realization
                if m == 0 && d>0 && forb_shocks[i_t]>prob_forb
                    forb = 1 #in forebearance!
                end
                i_f = forb + 1

                ###get consumption, wages, taxes, etc
                incwage = exp(γ_grid[i_e, 1] + γ_grid[i_e, 2] * θ + γ_grid[i_e, 3] * x + γ_grid[i_e, 4]*x^2) 
                w = incwage * (m == 1)
                taxes = tax(incwage) * (m == 1)
            
                c = interp_pol_nilf(θ_ind, i_e, i_x, i_a, i_t, i_d, i_f)
                if l == 1 && m == 0
                    c = interp_pol_unemp(θ_ind, i_e, i_x, i_a, i_t, i_d, i_f)
                elseif m == 1
                    c = interp_pol_emp(θ_ind, i_e, i_x, i_a, i_t, i_d)
                end

                ###determine payment made on loans
                periods_left = 10 - (i_t - 1) + 1
                periods_left += 1*(i_e==2) #update if somecoll
                periods_left += 3*(i_e==3) #update if college
                
                #bound
                periods_left = max(periods_left, 1)
                periods_left = min(periods_left, 10)
                
                #compute standard payment and bound
                pay_stand = payment(d, r_d, periods_left)
                pay_stand = max(pay_stand, p_min) #bound below
                pay_stand = min(pay_stand, p_max) #bound above by max payment
                pay_stand = min(pay_stand, d * (1+r_d)) #bound above by debt remaining
                debt_payment = pay_stand #initialize debt payment as the baseline "good standing" case

                #now check for delinquincy
                delinq = 0
                if (1+r)*a - pay_stand<0 && forb == 0 && m == 0 #not working, not in forebearance, and can't pay off loans with assets alone
                    delinq = 1
                    debt_payment = max((1+r)*a, 0)
                end

                #so now we have the baseline "good standing" payment. If individual is in forebearance, this goes to zero
                if forb == 1 || d == 0.0
                    delinq = 0 #not delinquent if in forebearance
                    debt_payment = 0.0
                end

                ###compute utility and record line of data    
                #utility
                util_leisure = λ_0 + λ_1*θ + λ_2*t + λ_3*(t-50)*(t>50) + λ_4 * t^3
                util = utility(max(c, cfloor), η) + util_leisure * (l==0)
                switch = 0

                #apply switching cost and update value of p
                if p==0 && l==1 && i_t!=t_start
                    util -= ξ
                    switch = 1
                elseif p>0 && l==0 && i_t!=t_start
                    util -= ξ
                    switch = 1
                end

                #NO LONGER convert taxes to NPV value
                line_simul = [θ, e, x, a, t, p, I, H, τ, s_choice, s_choice_val, prob_diff_cfact_I, prob_diff_cfact_W, l, m, w, taxes, c, util, switch, gov_exp, i_sim, i_v, d, debt_payment]
                line_num = (i_v-1) * nsim * nt + (i_sim-1) * nt + i_t #figure out what line of data we're in
                data_simul[line_num,:] = line_simul #update line

                ###evolution of state space. θ, e don't change. i_t will increment by itself.###      
                #increment experience
                i_x = i_x + 1*(m==1)
                
                #assets evolution
                ap = w - taxes + (1+r) * a - c - debt_payment
                ap = max(-cfloor, ap) #bound, just to be safe
                i_a = get_index(ap, a_grid)
                
                #debt evolution
                dp = d * (1+r_d) - debt_payment*(1-forb) #evolution of debt in good-standing or forebearance case
                if delinq == 1
                    dp = (d*(1 + r_d) - debt_payment) + 0.185 * (pay_stand - debt_payment) #add penalty to debt evolution
                end
                i_d = get_index(dp, d_grid)
    
                #update p AFTER having written the line
                p = 0
                if l == 1 && m == 0
                    p = 1
                elseif l == 1 && m == 1
                    p = 2
                end
            end
        end
    end
    data_simul #return simulated data
end
